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MATLAB 
تعليہ الماتلاب خطوة بخط وة‎ 


اعدا وتقدیم 


المهندس:- احمد محمد الفلاح الرابطى. 

التخصص:- الهندسة الكهربية والإلكترونية. 

القسم:- التحكم الألى. 

الجامعة:- (جامعة الجبل الغربي, جامعة الفاتح, جامعة الله آباد) 

المهنة والخبرة:- (معيد في كلية الهندسة, مهذدس في شركة ليبياذا, مهذندس في مصنع الأدوية 
بالرابطة, معلم ذي المعهد الصناعي بالرابطة, مهذدس منفذ ومتابع صيانة منظومة الدراسة 
والامتحانات في شعبية الجبل الغربي) 

Ahmad_Engineer21 @yaho0o.com -:JنigرiكJJا البريد‎ 

الم ستوى التعليم ي:- بك الوريوس ذي الهندسة الكهربي ة شعبة التحكم الألى من جامعة الجڊل 
الغربي. ودبلوما في الدراسات العليا في الهندسة الكهربية شعبة التحكم الألى من جامعة الذاتح. 
والآن دراسة الماجستير في جامعة الله أباد في الهند وتحضير لمناقشة رسالة الماجستير. 

السنة الدراسية:- 2010 - 2011 ف. 

الهواية:- المطالعة والشطرنج وكرة القدم. 

العنوان:- ليبيا - غريان - الرابطة. 
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24 99 


aes.) | emma | 


بسم الله الدى لا يحمد ولا يغذر ولا يسأل إلا هو وحده لا شريك له نعبده ولا ذشرك ډه 
شيء و صلى اللهم وسلم على النبي المصطفى خاتم النبدين, وعلى أله الطاهرين البررة وعلى 
أصحابه الأكرمين الدين نشروا الدين في البلدان وحملوا القران وحفظوا السنة, وعلى زوجاته 
الطاهرات أمهات المؤمنين وبعد... 

فمن خلال دراستي في لغة الماتلاب لاحظت أن هذاك عدة كتب تشرح البرمجة بلغة 
الماتلاب ولكن توجد ندرة في الأمثلة العملية في هده الكتب فنلاحظ الكاتب يكتذي بكتابة مال ام 
مثالين بسيطين قد لا يعطى طالب العم مراده وأيضا من خلال ملاحظتي لتخبط الكبير ابعض 
الط لاب في كتا ة الب رامج بلغ ة الم اتلاب ذي معمل الحاسوب. وعدم فهم كيف يتم ت صحيح 
الأخطاء . لهدا كتبت لكم قدر كبير من الأمظة العملية مع الخرج لتوصيل الفكرة ڊسهولة ويسر 
وبسرعة وبدون تعقيد. وقد تأكدت من النتائج للبرامج كلها في الحاسوب. . وكل هدا في سبيل 
تيسير العم فذسأل الله أن يجزيذا عن هدا العمل كامل الجزاء في يوم تزل فيه الأقدام اذه نعم 
المولى ونعم النصير. ونسأل كل من استفاد من هدا العمل الدى أخد منى ساعات طوال لتحضيره 
وإخراجه لكم على مثل هده الصورة المنظمة والواضحة أن يدعوا لذا في ظهر الغيب وذسأل الله 
القبول وعدم الرياء والنفاق فهو نعم المولى ونعم النصير . 


مفخرة للإنسان العلم 


واحذر يفوتك فخر ذاك المغفرس he‏ مرس كل فر قات 
من همته في مطعم أو ملببس ............واعلم بأن العلم ليس يناله 
في حالتيه عارياً أو مكتسي eee‏ | و ا ی و 
وافخو له ظببت الرقاة وكيس Sb Bes wll SS eesaniimo‏ 
RUE e SG I as E‏ 
اللذة في طلب العلم 
من وصل غانية وطيب عناف i abas‏ سهري لتنقيح العلوم ألذ لي 
احلى من الدوكاء والعشاف ١................9صریر‏ أقلامي على صفحاتها 
تفرك الفى الهل عن أواتى E E O O‏ 
في الدرس أشهى من مدامة ساقي a CIB lG es‏ 
نوما وتبغي بعد ذلك لحاقي انتا س قران ال دجا وتييتة 
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1- الجمع فى الماتلاب 


clear 


2- الطرح فى الماتلاب 


3- الضرب فى الماتلاب 


4- القسمة فى الماتلاب 
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elle 

clear 

a=4; 

b=5; 

c=; 

d=deconv (a,b) %SOr d=deconv (deconv (a,b) ,c) 
£=deconv (d, Cc) 


5-تمثيل الجدر والدالة الأسية واللوغاريتم الطبيعى والدوال المثلثية فى الماتلاب 


( (5*log10 (x) +2*x^2Z*sin (x) +sqrt (x) *lin (>) ) 


x=1; 
L=UECONT (ST LOGLU(R ETR ISIN FESOET (EI LOGS I ORO OTE I EIT AS 
LALOR) 


x=1; 
T= (SF LOGIU(E EDTA OTERO SCTE (RYT IOS Y7 (OSD (ORS ITE ASTD CLOG 
(x))) 


6- التفاضل فى الماتلاب 


Syms x 
f= ((xX5) (SX 4) H(A X3) = (2*72) = (GK) +9) 
O=OTEEFCE SX) 


Syms x 
LSC (ETI ESTE EAT EATS = (TE T= (O FE) F9) 
d=OTEF(EGA) 


Syms x 
£= (17 (1+%* 2) ) 
ASOTEFEF (EX) 


7 التگامل فی الماتلاب 


Syms x 
LZ (Ly (LEE 2)) 
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كل الخدت ااا 


2 ت 


٣ 


Syms x 
LS (RT IF ESR UAT EARS = ART T= (GFK ILO) 
d=1int (f, xX) 


Syms x 
T= (EDI E(SPL AT EAN I = (AK = (OTK) #9) 
d=int (f,1,2) 


SyYmMS xX a Bb 
f= ((x^5)+ (5*x^ 4) + (4*x^3) - (2*x^ 2) - (8 * x) +9) 
d=int (f,a,b) 


Example 1 


Fınd the total solution of the ODE 


ا 
df‏ 


subject to the ınıtıal condıtons KO) = 1 and ¥(0}) = - 


Solution: 


Yl ( > ke” E ke" 


(We must remember that the constants K, and K, must be evaluated from the total response). 


م. احمد محمد الفلاح الرابطى 


lo fınd the forced response, we assume a solution of the form 


38“ 


yr = de” 


{Ae "RAE "+ 3A" = 


Ahmad_engineer21 @yahoo.com 
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from which 4A = —-3 and the total solution 1s 
- 1 ب‎ d۴ 
YD) = yy+YyF = ke’ + ke 3e 


The constants F, and K; are evaluated from the given initial conditions. For this example, 


y(0) = 1 = Fe + Eye - 3e 


Î 
E TT 
dtl, _ , 
I 
Kk, - 3K; = 7 


vıelds KK = 2.3 and Kk; = LS. 


ا 


WF) = yy typ = 2.58 ° + L.3eê "3e 


Syms x L 

y=dsolve ('D2y+4*Dy+3*y=3*exp (-2*t)', 'y(0)=1', 'Dy (0)=-11') 
ezplot (yy, [O 4|) 

SPEEA 


9- حل معادلتين وثلاث معادلات للاآإخراج الثوابت 


Example l1 Example2 
<A | 
9k3+25k4=0 
) 0.68311۲3. r3=3 
25k3-19k44 0.393:1+3.81713=-1 
solution 

۲۲ 
Syms k3 k4 


E£I=T19*kKSFASTKA; 

LS NO E N 

[k3 k4|=solve(fl, f2 ) 

ells 

SVMS El TE FC 

EIS IEA 

a SOS CG E SG DERA 
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£3=0.393*r1 + 3.81/*r3+1 
[rl r2 r3|=solve(f1l, £2, £3) 


Exaımple3 


x" 2*y ^2+z=0 
x-(y/2)-alphatz=0 
xtzty=Û 


Solution 
ELE 
SVMS: KE. YF 3 AIPA 
2 0 
x¬ (y/2) -alpha+z=0 
X+Z2+Y=0 
[XVII] =SO1VE (REV LEE", TEVE I OID NAGE", GEST) 


0- إيجاد الجذور من معادلة متعددة الحدود والعكس 


Example 1 


Find the roots of the polynomial 


. 


Pp) = × 10 + 357 - 50x + 24 
P2) x) 


Solution: 


clc 

p1=[1 —-10 35 -50 24]; 
f£1=roots (P1) 

p2=[1 —-7 0 16 25 52]; 
£Z2=roots (P2) 


5 س کے سے سے سے کے س ج سے کے م سے س سے س ج س کے سے سے س کے س سے س کے کک سے کے س کے س‎ e e 


Example ا‎ 


1- It is known that the roots of a polynomial are J. 2. 3. and 4 Compute the coefficients of this poly- 


normal. 


2- It is known that the roots of a polynomial are —-J, —2, -3, 4+j5 and 4-j3. Find the coefficients 


of this polynomial. 


Solution: 


clc 

r1=[1 2 3 4]; 

£1=poly (£1) 

£2=[-1 -2 -3 -4+5[] -4-57 [; 
Ê£2=poly (F2) 


5 س سک س س س کے ب م س کس س س س سے سے کے س کے سے سے کے س کے س‎ e ب م س‎ e e mm em n 
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Example 3 


Evaluate the polynomial 


Clê 
p1=[1 -3 0 5 -4 3 2]; 
f£1=polyval (p1, -3) 


Example 4 


Let 


and 
6 4 2 
P= 2x x +4x + Ox l2 


1- Compute the product p, ‘ p; using the conv(a,b) function. 
2- Compute the product Pp,‘ Pp; using the [q,F] = deconv(c,d) function. 


Solution 


Cle 

p1=[1 0 -3 0 5 7 9]; 
p2=[2 -8 0 0 4 10 12]; 
f£1=conv (p1, P2) 

[q, r]=deconv (p1, P2) 


ج س س س س ت ت ت ن ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت کت 0 


Example 5 


Let “ 
2 - 8× + چ4‎ + [0× + 2 


ş usıng the PpOlyder(p) function.‏ مگ 


1- Lompute the derıvaive 
ر‎ Lompute the ıutesration ۳5 using the polyint(p) function. 


Solution: 
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clear 

p5=[2 0 -8 0 4 10 12]; 
f£1=polyder (p5) 
£2=polyint (p5) 


سس سس سے سے سے سے س سے سے سے سے کے سے سے سے کے سے سے سے کے کے کے سے کے سے سے کے سے سے سے سے سے سے سس 5 


Example 6 


R(x) mt _ (x - 4.8372x + 6.9971) + 0.6740x + 1.1058) + 1.1633) 
` Pgden (x -3.3520x + 3.0512)(x + 0.4216x + 1.0186)(x + 1.0000)(x + 1.9304) 


Pnum 
Fmd ~ ın polynomial form using the collect(s) function that 1s used to multiply two or more 
den 
symbolic expressions to obtain the result ın polynomual form. We must remember that the conv(pP:q) 


function 1s used wıth numeric expressions only, that 1s, polynomual coefficients. 


solution 


syms x 
pnum=collect ( (x^*2-4.8372*x+6. 9971) * (x^*2+0.6740*x+1.1058)( * (x+1.1633) ) 
pden=collect ( (x^*23.3520*x+3.0512) * (x^*2+0.4216*x+1.0186) * (x+1. 0000) * (x 
+1. 9304( )( 

R=pnum/pden 

pretty (R) 


سے سے سے سے سے سے سے سے سے سے سے س سے سے کے کے سے سے سے سے سے سے سے سے سے کے کے کے کے کے سے کے کے سس 5 


Example 7 


finds the residues, poles and direct term of a partial fraction expansion of the ratio of 


| 
| two polynomials B(s)/A(s) .If there are no multiple roots, 
| B()  R() RO) R(n) | 
س = س‎ 4+ -------- + „a. + ==------ + K(s) 
A(s) s-P(1) s-P(2) s - P(n) 
[R,P,K| = residue (B,A) 


bb ___ XxX^4+2xل^‎ 3-4 2+5+1 
ada xX^“5+4x4-2x*3+6x 2+ 2x+1 


| solution 


clc 

b=[1 2 -4 5 1]; 

a=[1 4 -2 6 2 1]; 
[R,P,K] = residue (b,a) 
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R= 0.2873 , -0.0973 + 0.17671, -0.0973 - 0.17671, 0.4536 + 0.00221, 0.4536 - 0.0021 
P = -4.6832, 0.5276 + 1.07991, 0.5276 - 1.07991, -0.1860 + 0.33651, -0.1860 - 51 


K =0‏ 
1- كيفية رسم الدوال فى الماتلاب 
Example .1‏ 


Write the MATLAB code that produces a simple plot for the waveform defined as 


7 


: 4 ت‎ RF | 4 
yY = Jf) = 3e  cosSt-2e sind 
7 


إ + 


ın the O ZF= 3 seconds interval. 


Solution: 


L=Û: UeOlLs S5 % Define t-axis in 0.01 increments 

VSR EDA Eya OOS FT EFT ODL FT Eo ENA TCS E TC 
(t+1) 

pot (EY; 

GET, 


xlabel ('t'); 
ylabel ('y=f(t)'); 
LLLLS (PLOT TOF Examogle A.131") 


Example 2 


Plot the functions 


3 


2 2 ENE 2 
J = FIMX, Z = COS XX, W = SIM X* COSY, V = SIR XxX COS" 


ر 
3 


in the interval O0 = xX & 21 using 100 data points. 


1-Use the plot command to display these functions on same windows on the same graph. 


2- Use the SUbplot command to dısplay these functions on four windows on the same graph. 


Solution: 


1-Use the plot command to display these functions on same windows on the same graph. 


Cle 
clear 
x=linspace (0, 2*pi, 100); % Interval with 100 data points 


y=(sin(x).^ 2); 
z= (cos(x).^ 2); 
W=y.* Z; 


v=y./ (z+teps); % add eps to avoid division by zero 
plot (x,y,'b',x,Z,'g',x,W,'E',x,v,'y'); 
grid on 


axis ([0 10 0 5]); 


ی کک م ی سے سے ج ج ص ب س م س س ج ج س س ج ب س ب e‏ س س س ج ب ص ج ص س س دک e e‏ د n‏ س ج س س e‏ و و ک س م س م س م س س س ج س س ج ج ج ب کے س 5 
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2- Use the SUbplot command to display these functions on four windows on the same graph. 


Clc 
clear 
X=1INSDSCE(O; 2 ODI, 100); < Interval with 100 data points 


SS CR ET 
Z=(COS (XxX). 2); 
W=Y.* ZF 

V=y./ (Zz+teps); 
SUDPLOL (ASIL 
PLO (Xr Y; 

AXKISCLVU Z5I U LIF 
title ('y=(sinx) ^2' 
GELA OF 

SUDPLOT SLA), 
BLOC A2; 

ARLISTLU ZBL O Ll) 
LILLE E= (COK AT 
GEL O0 

SUDPBLOL LES, 

plot (x,w); 

axi1iS5([]Û0 2*51 O0 0.31); 

title ('w=(Ssinx) *2* (COSX)  2'); 


o0 


UGE LETE OF TOUE SUDPBLOLS 


o0 


Upper FIORE GF TOU SUDDIGES 


7 


o0 


lower left of four subplots 


GEL. OM 
SUBPLOT, > LOWEE FIGRL OF TOUE SUDDLOLS 
PLIOEIRE TI, 


ax1S5([0 2*5i 0 4001; 
LILLE TS OINE LT OO TDS 
GEIG BR 


E 2=(£05#) 
1 1 
1.5 5 1 
1 1 
1 2 1 2 4 B 
w= (sinc f sincosx)“ 
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مثال :- اكتب برنامج يحسب جدول الضرب ويعرضه فى شكل منظم باستعمال الأمر إه؟؟ 


for j=1:10; 


b=b+1; 
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تمتي ل المصفوفات فى المات لاب 
1- العمليات الحسابية للمصفوفات فى الماتلاب() ((Matrix Operations‏ 


Example L.1 


Compute A + B and A - B given that 


3 and B = ك‎ 
4 25 


Solution: 


AL 2A SF O1. AI < Define matrices A 
B=] 3g واو 2 اك .و0‎ < Define matrices B 
m1=A+B <> Add A and B 
m2=A-B < Subtract B from A 


Example C.2 


Multıply the matrıx 


Solution: 


(- 3+ j2) x 1 ا‎ 


+74 9+76 


-3+72(×»2 )-3+[2( ×3 
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Check with MATLAB: 
LE 
clear 
k1=5; <> Define scalars kl 
k2= (=3 ± 2(7; < Define scalars k2 
وا 2 وض ل‎ <> Define matrix A 
ml=k1*A MULLIS MALCEIS A By Conc lant EI 
mM42=K2*A SMULCIDILIYS MAEEI A GS CONSTR KL 


Example C.3 
Matrices C and ÛD are defıned as 


1 
> 2ا‎ 3 4 and D = |_7 


e 


اا 


Compute the products C'۰. D and D.C 


Solution: 


The dimensions of matrices C and D are respectively I x 3 3 x I ; therefore the product C. D ıs 


feasıble, and wıll result ın a J x J , that ıs, 


CS Cal O UCN CD O a 


“j 


ر 


The dımensions for D and C are respectively 3 x I I1 x 3 and therefore, the product D - C 1s also 


feasible. Multiplication of these wıll produce a 3 x 3 matııx as follows: 


(1)‘(2) (D(3) (G)°(4) 
(1)‘(2) (1° (3) (1)°(4) 
)2()2( )2()3( )2()4( 


C=[2 3 4|; <> Define matrices C and D 
D=[1; =1; 2|; %4 Define matrices C and D 
ml1=C*D eS MuOlTIGIT CEC BT DB 
M2 =D*C E MulLEILBLF D BGC 
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((Determinants of Matrices )) تlİوêصملل حساب المحددات‎ -2 


Example L.4 


Matrices A4 and BÊ are defıned as 


and FB = ٌ 
: د‎ 


CLompute detA and d{etB. 


Solution: 


detlAd 
detB 


B= =F 2 O1; <> Define matrices A and B 
det (A) %< Compute the determinant of A 
det (B) < Compute the determinant of B 


Example L.5 


Lompute detA and defB ıf matnces 4A and B are defined as 


7 3 


2 لا ل 
0 3 0 


Solution: 


det(A)= 2(0*0-1*1) - 3(1*0-1*2) + 5(1*1-0*2) = -2+6+5 = 9 


det(B) = 2(0*(-6)-(-2)*(-5)) - (-3)(1*(-6)-0*(-2)) +4(1*(-5)-0*0)= -18 


Check with MATLAB: 


A= A2. 3 SE 1L OL; 2L UI <> Define matrix A 
=| وغ وت ل رت ل1 لھ وھ‎ < Define matrix B 


det (A) < Compute the determinant of A 
det (B) < Compute the determinant of B 
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((Cramer’s Rule )) قاع‎ -3 


Let us consider the systems of the three equatıons below 
Appx F ApaY TF Aj3 ار کک‎ 
FTES + FE + شوو‎ 
and let 
ی‎ A7] 
Dı =| B ao, 
C a3 


a. 
A 


x 


provıded that the determınant A (delta) ıs not zero. 


Example C.10 
Use Cramer’ rule to fınd vy, Yڊ,‎ and vy ıf 
و3۷ + 3-۷ - ر2۷‎ = 0 
گ = ر۷ ر۷ و۷‎ 
ر۷‎ + 3۷-4-۷ = Û0 
and verıfy your answers with MATLAB. 
Solution: 
Rearrangıng the unknowns Y and transferring known values to the right side, we get 
و = و۷ک + ر۷ ر۷‎ 


کے کک و2۷ ج 3۷7 ج HV;‏ 


Now, by Cramer's rule, 


|3 ي : 7 2# ا4 ل ت 
_g = ]5+8+24+ 30+ [10-8 = 835‏ ۾„ |37 3_ = A =| _y _3 7| 4 73 = F+06-[]2+27+4+4 = 35 D0,‏ 
ا ا کک کک 1 ت ؟ 3 : 
۶7 کې اګ 7ر ي 7[ 3 7 1[ 3 
ر 3 ل 2 
70- = 16-30-48-72+16-20- = £ 4 |27 £ 4= 
4 3 أ1 4 31 
4d‏ ل ړ 4d‏ 
55 = 0[-ف[- 43 +24-24-20- = چ ي إميم 3 ب4_ |= 
1 3 4 1 3 
Then,‏ 
TT Gg aT‏ 
A 33‏ 7 35 ۸ ت 7 35 A‏ | 
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We will verify with MATLAB as follows. 


clear 
% The following code will compute and display the values of v1, v2 
and v3. 


o9 


B=[2 -1 3;—-4 —-3 -2; 3 1 -1[; The elements of the determinant 
D of matrix B 
delta=det (B) ; 


matrix B 


d9 


Compute the determinant D of 


d1=[5 —-1 3; 8 -3 -2; 4 1 -1[; % The elements of D1 

detd1=det (d1) ; % Compute the determinant of D1 
d2=[]2 5 3; -4 8 -2; 3 4 -1[; % The elements of D2 

detd2=det (d2) ; % Compute the determinant of D2 
d3=[2 -1 5; -4 -3 8; 3 1 ;إ4‎ % The elements of D3 

detd3=det (d3) ; % Compute he determinant of D3 
v1=detd1/delta; % Compute the value of v1 
v2=detd2/delta; % Compute the value of v2 
v3=detd3/delta; % Compute the value of v3 

disp ('v1=');disp (v1); % Display the value of v1 

disp ('v2='); disp (v2); % Display the value of v2 

disp ('v3=') ; disp (v3) ; % Display the value of v3 


5 س کے س ج بس سے ب کے ج کے ج ج س سے س ج س س ب ج ب مح س ص س س س ی س ب م ب ج س س م س س سے س کے س کے سے سے سے کے س کے س کے س سے کے س کے س‎ e e 


((The Inverse of a Matrix )) dوفصملا 4-حساب معكوس‎ 


Example C.14 


Matrıxz A ıs defined as 


Lompute ıts inverse, that ıs, find 4 


Solution: 
Here, detA = O+8S+ 12—-9-JIO—-û = —2, and sınce thıs 1s a non-zero value, ıt ıs possible to com- 


pute the inverse of 4 using 


Then, 
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Check with MATLAB: 


clear 

A=[1 2 3; 1 3 4; 1 4 3]; 
invA=inv (A) 

%format long; invA 
%format short; invA 


5 ب م ی س س س کے س‎ e e ب ص ا ص س و‎ e ب‎ a e e n a n e ا‎ n ی‎ e e: e e e e د د کے ب‎ e س و‎ e ص ج ج ب س ب‎ e e e e e n e e: a r 
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(( Solution of Simultaneous Equations with Matrices )) ٽIİوذصملا حلول المعادلات الأنية باستخدام‎ -5 ) 


Example L.16 


For the system of the equations 


Ll - 97 
XxX +t 2x” + 3X; = 
|×, ++ و2‎ = 8| 


compute the unknowns XJ, XxX, and Xx; usıng the ınverse matrıx method. 
Solution: 


În matrıx form, the gıven set of equations 1s AX = B where 


9 ر 7 3 2 
ما =8 ,را =& ,إو 7 [اه4 
- 4 
X3‏ ر آ7 3 
Then,‏ 
X= AB‏ 
LENE‏ 
و 7 3 2 
1 کے ا يات 
ی ر لر هھ 


و ا ا 

detdA = IS8 and adjA= | 7 7J 5 

ر ”ر ا 

Therefore, 

د 

d ف‎ add —= ا‎ 5 

3 1 7 ور > ل وی = ^ 

ر و اک 

we obtaın the soluton as follows. 

1 3 79 33 357 18 4 
E N N A SS 
8ا" ” ک5‎ 5 37 18 0.28 
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Check with MATLAB: 


clear 
A=[2 3 1; 1 2 3; 3 1 2]; 
B=[9 6 8]'; 

X=A\B 

M=inv (A) *B 


س 
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Example C.17 
For the electric cırcuut of Figure C.1, 


£2 2 2 ل 


و3 ا ا 17 ا 
a 3‏ 
1 ا FOOTY ٤ ١‏ = ¥ 
I?‏ > | 9 3 و : 


Figure C.I. Circuit for Example C.l7 


the loop equations ate 


IOF, 9F, = 100 
97, + 201-91; = 0 
9I + 15I = 0 


Use the ınverse matrix method to compute the values of the currents I, , I , and I, 


Solution 


For thıs example, the matrıx equatıon ıs RI = FP or I = R7 F, where 
10-9 0 100 1; 

R= |_g 20 gj, F=|j pj and I=, 

0 75 0 Ff 


E eT 8 
The next step 1s to find RF . Thus ıs found from the relation 


-[ 1 


` detR 


Iherefore, we fınd the determınant and the adjoınt of R . For thıs example, we fınd that 


f9 735 87 


detR = 975, adjR = 135 150 00 
81 90 119 
Then, 
219 135 81 
E 
8I1 00 119 
and 
1 1 22.46 
1= إ0 |90 150 133| = ر1‎ ^ g7 133 > |35 
1; S1 00 110| 0 81 8.31 
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Check with MATLAB: 


clear 

R=[10 —-9 0; —9 20 —-9; 0 —-9 15]; 
V=[100 0 O]'; 
I=R\V; 

disp ('I1='); 
disp (1 (1)) 
disp ('12='); 
disp (1 (2)) 
disp ('I13='); 
disp (1 (3)) 
M=inv (R) *V 
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Example L.18 


For the phasor cırcuut of Figure C.18 


Figure C.3. Circutlt for Example C.1838 


the current Jy can be found from the relation 


F—Fڊ‎ 

څک ڪ= م[ 
R‏ 4 
3 

and the voltages F, and F, can be computed from the nodal equations 


F, IO FF, F0 
پک‎ f 
85 100 ` j200 


ancl 


F۳ E 170° ۳” ۸ F۳” 2 ل‎ 
س‎ e ا‎ 
[100 100 1 


Compute, and express the current I, ın both rectangular and polar forms by first sımplifyıng lke 


terms, collectıng, and then wrıtıng the above relations ın matrix form as FF = I, where 


Y = Admittance, WV = Voltage , and I = Current 


Solution: 


The YT matrix elements are the coefficients of 7 and Fr. sımplifyıng ancl rearrangıng the nodal equa- 


UONS , we get 


(0.0218 - 70.003) PF, 0.01¥; = 2 
—O.OIP, +(0.03 +7O0.OI)F; = JL7 


Next, we write ın matrix form as 


0.0218 -j0.005 01 7ı 2 
-0.01 0.03+j0.01| |, = |7 


4 F7 I 
where the matrices F, F, and J are as ındıcated. 
V1 = 1.0480e+002 + A4.8448e+0011 


Vê = 53.4162 + 55.34391 


Iherefore, ın polar form 


I, = 0.518 2-6.33° 
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Check with MATLAB: 
Cle 
clear 
Y=[0.0218-0.005] -0.01;-0.01 0.030.017 [; % Define Y, 
I=[2; 1.73]; % Define I, 
VY T} % Find V 
M=inv (Y) *I; 
fprintf ('\n'); % Insert a line 
disp('V1 = '); 
disp (V(1)); % Display values of V1 
disp ('V2 = '); 
disp (V(2)); % Display values of V2 
R3=100; 
IX= (V (1) -V (2) ) /R3 % Compute the value of IX 
magIX=abs (IX) % Compute the magnitude 
of IX 
thetalX=angle (IX) *180/pi % Compute angle theta in 
degrees 
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Example I1. 


Sımplify the complex number z and express ıt both ın rectangular and polar 
Torin. 


)3 + 76 ()1 + 72 ( 


SoluifiOon: 


% Evaluation of Z7 


% the complex numbers are entered 


N 
ى‎ 
1 


3+4*]; 
22 = 5+2*[; 
theta = (60/180) *pi; 
Z23 = 2*exp (j]J*theta) ; 


d9 


angle in radians 


Z24 = 3+6); 
Z25 = 1+2*j; 

% Z_ rect is complex number Z in rectangular form 
disp('Z in rectangular form is'); % displays text inside brackets 


Z rect = Z1*Z2*Z3/ (Z4+Z5) 


Z_mag = abs (Z_rect); % magnitude of Z 
Z_angle = angle (Z_rect) * (180/pi) ; % Angle in degrees 


disp ('complex number Z in polar form, mag, phase'); % displays text 
%inside brackets 

Z_polar = [Z_mag, Z7_anglej] 

diary 
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Example 1. 


Write a function file to solve the equivalent resistance of series connected re- 


sistors, R1, R2. R3, .... Rn. 


Solution: 


E EE EE E E EEE EEE E EE RE E E EEE ERE E E EE ER E EEE E EEE 
function req = equiv_sr (EF) 
% equiv _ sr is a function program for obtaining 
% the equivalent resistance of series connected resistors 
% usage: req = equiv _ sr (F) 
% Fr is an input vector of length n 


% req is an output, the equivalent resistance (scalar) 


n = length (r); % number of resistors 
req = sum (FEF); % sum up all resistors 
end 


The above MATLAB script can be found in the function file 
equiv_sr.m, which 1s available on the disk that accompanies this book. 

Suppose we want to find the equivalent resistance of the series 
connected resistors 10, 20, 15, 16 and 5 ohms. The following statements 
can be typed in the MATLAB command window to reference the 
function eqU1IVv_SF 


a = [10 20 15 16 5]; 
Rseries = equiv_sF (a) 


diary 


The result obtained from MATLAB 1s 


Rserles = 
66 
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Example 1.. 


Write a MATLAB function to obtain the roots of the quadratic equation 


ax + bx + Cc = 0 


YOIUTION: 


ج س س ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت کت کک کک ت 0 


function rt = rt_ quad (coef) 


% rt_ quad is a function for obtaining the roots of 


% of a quadratic equation 
% usage: rt = rt_ quad (coef) 


% coef is the coefficients a,b,c of the quadratic 


% equation ax*x + bx + Cc =0 


% rt are the roots, vector of length 2 


% coefficient a, b, c are obtained from vector coef 


a = coef (1); b = coef (2); c = coef (3); 


int = b^2 —- 4xax*c; 
if int > 0 
srint = sqrt (int); 
x1= (—-b + srint) / (2*a); 
x2= (-b - srint) / (2*a); 
elseif int == 
x1= —-b/ (2*a); 
X2= x1; 
elseif int < 0 
srint = sqrt (-int); 
p1 = -b/ (2*a); 
P2 = srint/ (2*a); 
x1 = pl+tp2*]; 
x2 = p1-p2*j; 
end 
rt =[x1;x2]; 
end 


We can use m-file function, rt_quad, to find the roots of the following 


quadratic equations: 


(a) x +3x+2=0 (bJ) £ +2x+1 =0 (¢) £ -2x +13 =0 


م. احمد محمد الفلاح الرابطى 


%Sdiary ex1.dat 
ca = [1 3 2]; 


ra = rt_ quad (ca) 
cb = [1 2 1]; 

rb = rt_ quad (cb) 
cc = [1 -2 3]; 
rC = rt_ quad (cc) 
diary 


Ahmad_engineer21 @yahoo.com 
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aaa DS 


clc 
close all 


a = input( ' a = '); 
b = input( ' b = '); 
c = input( ' cC = '); 


x1 = ( - þ + sqrt( b^2-4*a*c) ) / (2*a) 
x2 = ( - þ + sqrt( b^2-4*a*c) ) / (2*a) 
1f imag (x1)==0 & imag (x2) == 


1f x1==x2 
str=' ident ' 
else 
str= 'real' 
end 


elseif real (x1)==0 & real (x2) == 
str=' imag ' 
else 
str=' comp ' 
end 
bigstr=[' (x1=',num2str (x1), ')—--',' (x2=',num2str (x2), ')--' ,ster]; 
msgbox (bigstr) 


مثال 
برنامج لقياس الوقت الي تستغرقه للوصول لبلد علي بعد 800 كيلومتر يعني اننا سندخل 


طريقة المواصلات هل هى عربة أم حافلة أم طائرة ,العربة تسير بسرعة 120 
كيلومتر إساعة والحافلة بسرعة 860 كيلومتر /ساعة والطائرة بسرعة 200 كيلومتر إساعة 


clear 
CLC 
close all 
a=input ('enter your transportation method :','s'); 
switch a 
case 'car'" 
t=800/120 
msgbox (['your trip will take ',numZ2str (t),' hours']); 
case 'bus' 
t=800/80 
msgbox (['your trip will take ',numZ2str (t),' hours']); 
case 'plane' 
t=800/200 
msgbox (['your trip will take ',numZ2str (t),' hours']); 
otherwise 
msgbox ('inter valed tm') 
end 
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1۔ تمثیل (( sمpہا‏ اہ ٥ط‏ )) فی الماتلاب 


Repeating with for loops 
Syntax of the for loop 15 shown below 
for k = array 


commands 
end 


The commands between the for and end statements are executed for all values stored 1n the 
AarTAV. 


Example 1 


Suppose that one-need values of the sıne function at eleven evenly spaced poınts rn/10, for 
ı1 = 0, 1, ..., 10. To generate the numbers iıı question one can use the for loop 


Solution 


for n=0:10 
x(n+1) = sin (pi*n/10); 

end 
x 
clc 
H = zeros(5); 

for k=1:5 

for l=1:5 
H(k,1) = 1/ (k+1-1); 
end 

end 
H 
clc 
A = zeros(10); 
for k=1:10 

for l1l=1:10 
A(k,l1l) = sin (k) *cos (l1) ; 
end 

end 
k 1:10; 


sin (k) '*cos (k); 


سے سے سے کے سے سے سے سے سے سے سے س سے سے کے کے سے سے سے سے سے سے کے کے کے سے سے سے سے سے سے سے سے کے کے سے ست ست سے سے سے کس 5 
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2۔ تمثیل ((ئp٥٥1‏ ٥اذطس‏ ١ط)))‏ فی الماتلاب 
Repeating with while loops‏ 


Syntax of the while loop 1s 


while expression 
statements 
erd 


Thıs loop 1s used when the programmer does not know the number of repetitions a prior1. 


Example 1 
This process 1s continued till the current quotient 1s less than or equal 


| to 0.01. What 1s the largest quotient that 1s greater than 0.01? | 


| Solution | 


q = q/2; 
end 


3- تمثیل (( sرہناcںconstru ))the if-else-end‏ فی الماتلاب 


The If-else-end construction s5 


Syntax of the simplest form of the construction under discussion 1s 


if expression 
com mMmMaAnûÛs 
end 


This constructHon 1s used 1f there 1s one alternative only. Two alternatives require the construction 


if expression 

commands (evaluated if expression is true) 
else 

commands (evaluated if expression is false) 
erndl 


If there are several alternatives one should use the following construction 


if expressiaorıl 

commands (evaluated if expression 1 is true) 
elseif expression 2 

commands (evaluated if expression 2 is true) 
elseif ... 


else 
commands (executed if all previous expressions evaluate to false) 
emd 
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Example 1 


Chebyshev polynomials T(x), nı = 0, 1, ... of the first kınd are of great importance ın numerical 
analysis. They are defined recursively as follows 


Ta(x) = 2ZxT,-1(S) - T-2(D, =2, 3, ..., ToD =1, TS =x 


Implementation of thıs definition 15 easy 


Solution 


س س س ت س ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت کے سک کت کت کت کت کک کت کک کک ت 0 


function T = chebt (n) 
% Coefficients T of the nth Chebyshev polynomial of the first 


kind. 
% They are stored in the descending order of powers. 
tO = 1; 
t1 = [1 0]; 
if n == 
T = t0; 
elseif n == 1; 
T = t1; 
else 
for k=2:n 
T = [2*t1 0] —- [O0 0 t0]; 
tO = t1; 
EL ë5 TT; 
end 
end 
اسند عاء‎ 
ELE 
n=3 
coff = chebt (n) 
diary 


Thus T(x) = 4F -3x 
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4- تمثیل ))the switch-case construct: os((‏ فى الماتلاب 


The switch-case construction 


Syntax of the switch-case constructHon 15 


switch expression (scalar or string) 
case valuel (executes if expression evaluates to valuel) 
commands 
case valued (executes if expression evaluates to value2) 
commands 


otherwise 
statements 
erd 


Switch compares the input expression to each case value. Once the match 1s found 1t executes the 
associated commands. 


Exaınple 1 


In the following example a random integer number x from the set {1, 2, ... , 10} 1s generated. If 
x =1 or x =2, then the message Probability = 20% 1s displayed to the screen. If x = 3 or 4 or 3, 
then the message Probability = 30% 1s displayed, otherwise the message Probability = 50%0 1s 
generated. The scrıpt file fswitch utilizes a swıtch as a tool for handling all cases mentioned 
above 


Solution 


% Script file fswitch. 


x = ceil (10*rand) ; % Generate a random integer in {l1, 2, ... , 10} 
switch x 
case {1,2} 
disp ('Probability = 20%'); 


case {3,4,5} 


disp ('Probability = 30%'); 
otherwise 
disp ('Probability = 50%'); 


Note use of the curly braces{ }after the word case. This creates the 
so-called cell array rather than the one-dimensional array, which requires 


use of the square brackets|[]. 
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(( Rounding to integers. Function ceil, floor, fix and round )) پيرقتll‎ Jلاود‎ -5 


We have already used two MATLAB functions round and ceil to round real numbers to integers. 
They are briefly described ın the previous sections of this tutorial. A full lst of functions 
designed for rounding numbers 1s provided below 


۹C floor Round towards minus 1nfinity 
__ round 


Round towards nearest inteseTr 


Example 1 


To illustrate differences between these functions let us create first a two-dimensional array of 
random numbers that are normally distributed (mean = 0, variance = 1) using another MATLAB 
function randn 


Solution 


randn ('seed', 0) % This sets the seed of the random numbers 
% generator to zero 
T = randn (5) 
A = floor (T) 
B = ceil (T) 
C = fix(T) 
D = round (T) 
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Example | 


In the following m-file functions floor and ceil are used to obtain a certain representation of a 
nonnegative real number 


Solution 


function [m, r] = rep4 (x) 
% Given a nonnegative number x, function rep4 computes an integer m 
% and a real number r, where 0.25 <= r < 1, such that x = (4^m) *r. 
if x == 
m = 0; 
r= 0; 
return 
end 
u = 1log10 (x) /log10 (4) ; 


if u < 0 


m = floor (u) 
else 
m = ceil (u); 
end 
r = x/4^m; 
اسند عاء‎ 
clc 


[m, r] = rep4 (pi) 


ج س س ت ت ت ت س س ت ت ت ت ت ت ت ت ت ت ت ت ت ت ن ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت کت کت کے کت کت کے کت کے کے کت کے کے کت کت کے کے کت کے کت کک کک کک کت کک کے 0 
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((MATLAB graphics )) بںتlnll 1-الرسم فى‎ 


Example 1 


In this example the graph of the rational function f (x) = 
1+ 
using a variable ıuımber of points on the graph of f(x) 


Solution 


% Script file graph1l. 
% Graph of the rational function y = x/ (1+x^2) . 
for n=1:2:5 
n1O0 = 10*n; 
x = linspace(-2,2,n10); 
y = x./(1+x.^2); 
plot (x,y, 'r') 
title (sprintf ('Graph %g. Plot based upon n = %g points.', (n+1)/2, n10)) 
axis ([-2,2,-.8,.8[) 
xlabel ('x') 
ylabel ('y') 
grid 
pause (3) 


% Script file graph2. 
% Several plots of the rational function y = x/ (1+x^2) 
% in the same window. 
k = 0; 
for n=1:3:10 
n1O0 = 10*n; 


* 
1 


linspace (-2,2,n10) ; 
x./(1+x.^2); 
k = k+1; 


1 


subplot (2,2, k) 

plot (x,y, 'r') 

title (sprintf ('Graph %g. Plot based upon n = %g points.', k, n10)) 
xlabel ('x') 

ylabel ('y') 

axis ([-2,2,-.8,.8[) 

grid 

pause (3) ; 


end 
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Exaınple 2 


Using command plot you can display several curves 1n the same Figure Window. 


We wıll plot two ellipses 


x-7) (y-8 
= 4= ° =[ an ا ) ا‎ ) =1 
6 81 4 36 


using cominand plot 


Solution 


-2 + sin (Ct) 


X(T) = 3 + 6cCcOs(t), Y(E) 


1 
“J 


x(t) + ACOS(L), Y(L) = 8 + bslin(t). 


% Script file graph3. 
% Graphs of two ellipses 
% x(t) = 3 + 6cos(t), y(t) = —-2 + 9sin(t) 


% and 


d9 


7 + 2cos(t), y(t) = 8 + 6sin(t). 


x(t) 
t = 0:pi/100:2*pi; 

x1 = 3 + 6*cos(t); 

y1 = -2 + 9*sin(t); 

X2 = 7 + 2Z*cos(t); 

y2 = 8 + 6*sin(t); 

plot (xl1,y1l,'r',x2,y2,'b'); 

axis ([-10 15 -14 20[) 

xlabel ('x') 

ylabel ('y') 

title ('Graphs of (x-3)^2/36+(y+2)^2/81 = 1 and (x-7)^2/4+ (y-8)^2/36 =1.'( 


grid 
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vellow 
magenta 


creen 
blue 

white 
black 


Exaınple 3 


If function axis 1s not used, then the circular curves are not necessarily circular. To justify this let 
us plot a graph of the unit circle of radius 1 with center at the origin 


Solution 


t = 0:pi/100:2*pi; 
x = cos(t); 


sin (t) ; 


1 


plot (x,y) 


% Script file graph4. 


% Curve r(t) = < t*cos(t), t*sin(t), t >. 
t = -10*pi:pi/100:10*pi; 
xXx = t.*cos(t); 


y = t.*sin(t); 

plot3(x,y,t); 

title ('Curve u(t) = < t*cos(t), t*sin(t), E >') 
xlabel ('x') 

ylabel ('y') 

zlabel ('z') 

grid 
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مثال 
نقوم برسم وتحميل عدة رسمات فى نفس الشكل ونوضح كل خط بلون معين وقيمته على | 


clc 
close all 
x=linspace (0,10,1000); 
a=.1:.1:.6; 
c='b rt mMmMCEK YY"; 
for i=1:6 
y=sin (x) .*exp (-a (i) *x) ; 
plot (x,y, c(i) ) 
hold on 
end 


legend('a=.1','a=.2','a=.3','a=.4','a=.5','a=.6'( 


Figure 1 ١ E 3 


Fila Edit Yiew Insert Tools Desktop Windom Help ل"‎ 


DCEHã|R| AAC ® | Bj a A 
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Example 


Find first and second derivatives for F(<)=x<^2+2x<+2 
Solution 


To find first and second derivatives of Pn (x) -———- )‏ س 


clc 

a=[1 2 3]; 

syms x 

p=a (1); 

for i=1; 

p=a (i+1) +tx*p; 

end 

disp('First derivative') 
p2=p+x*diff (Pp) 

disp ('Second derivative') 
p22=diff (P2) 


First derivative 


p2 = 
2+2 ×× 


Second derivative 


P22 = 


م. احمد محمد الفلاح الرابطى Ahmad_engineer21 @yahoo.com‏ 


الى كل من استفاد من هدا الجهد لا نسألكم الشكر ولا الثناء انما نسالكم دعوة صادقة في جوف الليل عسا ان تنفعنا في يوم تزل فيه الأقدام 


Example 
P4(x<)=3x<^4-10x^3-48x^2-2x<+12 at r=6 deflate the polynom1al 
with Horners algorithm Find P3(X). 


| Solution 


HOENEE alOGOEFITI =a‏ ت ج ت ی 

clc 

a=[]3 -10 -48 -2 127; 

r=6; 

b(1)=a(1); 

p=0; 

n=length (a) ; 

for 1=2:n; 
b(1)=a(1)+r.*b(1-—1) ; 

end 

syms x 

for 1=1:n; 
p=p+b (i) *x^ (4-i) ; 

end 

disp ('P3 (x) =') 


3*£*3+8*x× 2-2 


Numerical Integration 
1- Trapezoidal Rule 


uu‏ ا 

1 X 

xX =a x=b * 1 n1 a1 Xn 
f 


Figure . T[rapezoidal rule. Figure. Composite trapezoidal rule. 


The composite trapezoidal rule. 
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H—~] 
I= > i= [f(x +2 f(2) + 2f) +۰. +2 f(1) + fm] 


1=1 


| Example 
Suppose we wished to Integrate the function trabulated the table 


below for Jf (x)=ge* over the ınterval from x=1.8 to x=3.4 using n=8 


b 3.4 
Am=[ f e ٢ (e* dx 


kL 6.050 | 7.389 | 9.025 | 11.023 | 4 20.086 29.964 | 36.598 | 4.1 
ا‎ 


Solution 


== TEaDOZOLISL Rl 
clc 
a=1.8; 
b=3.4; 
h=0.2; 
n= (b-a) /h 
£=0; 
x=2; 
for i=1:n; 

%c=a+t (1-1/2) *h; 

%Sf£=f+ (c^2+1) ; 

f= (f+exp (x) ) 

x=x+t+h; 
end 
Am_approx=h/2* (exp (a) +2*f+exp (b) ) 
syms +t 
Am exact=int (exp (t),1.8, 3.4) 
error=Am exact—Am approx 
E_t= (error/ (Am_approx+error) ) *100 
E_a= ( (Am _ approx—Am_ exact) /Am_ approx) *100 


5 س ج ب سے س وک سے سے س و س سے ص سے وک س کے ا‎ e a ا و‎ e e e o e a an e e 


2-Simpson’s 1/3 rule 


| 
X 8 


1 


b 


ٍ ا‎ Figure, Composite Simpson's 1/3 rule. 
Figure . Simpson's 1/3 rule. 
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The composite Simpson’s 1/3 rule 


bh ا‎ e 8 
IE AX x [I = [1 (x1 JF 4 | 2) ۳ ر‎ (Xal 4 f)x4( E 2 (Xn-a2) F 4 fl Xn-1 IF (Xn) 3 


f 


| Example 
Suppose we wished to integrate the function using Simpson’s 1/3 rule and 


Simpson’s 3/8 rule the table below for 1 (x)= ge* over the interval from x=1.8 


b 3.4 
Am=[ f e ا‎ (e dx 


to x=3.4 using n=8 


x» 16 18 2 22 2.4 | 2.6 28 | 3 32 34136 38 


$==-SINMpPSON SS 1/3 Ilo a 
clc 


Hh Š5 bD' Û © 
U IU IU I Il 
O~oOwWpH 

N #* o0 


Ba) 


e 


m=0; 
for x=2: (hth) :3.2; 
f= (f+exp (x) ) ; 
end 


for x=2.2: (hth) :3; 
m= (m+exp (x) ) ; 

end 

Am_approx=h/3* (exp (a) +4*f+2*m+exp (b) ) 
syms +t 

Am _ exact=int (exp (tt), 1.8, 3.4) 

error=Am exact—Am approx 

E_t= (error/ (Am_approx+error) ) *100 

E_a= ( (Am _ approx—-Am_ exact) /Am_approx) *100 


ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ت ت ت ت ت ت س 


3-Simpson’s 3/8 rule 


The composite Simpson’s 3/8 rule 


3 


j 
flxlax a 1= [flx) +3 faa) + 3 faa) + 2fl) + °°° +3 F(a) + 3 fn) + FAD) 
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a E GEG EEG Simpson's 3/8 rule —-——-—- 


h=0.2; 
n= (b-a) /h; 


for x=2:h:2+h; 
f=f+exp (x) 


x=x+t+h; 


m=exp (x) ; 


for x=2.6:h:2.6+th; 
f=f+exp (x) ; 


end 


x=x+h; 
m=m+texp (x) ; 
x=x+h; 


f=f+exp (x) ; 


syms tt 
Am exact=int (exp (tt), 1.8, 3.4) 
error=Am exact—Am approx 
E_t= (error/ (Am_ approx+error) ) *100 
E_a=( (Am approx—-Am_ exact) /Am_approx) *100 


کے کے کے کے کے سے سے سے کے کے سے سے سے کے سے سے سے کے سے کے سے کے سے سے سے سے کے کے کے سے سے سے کے سے کے سے کے سے کے سے کے کے کے کے سے کے کے کے کے کے کے کک 2 


© 


clc 


a=1.8;b=3.4;h=0.2; n= (b-a) /h; f=0;m=0; 


for x=2:h:3.2; 
switch x 
case {2,2.2} 
f£=f+exp (x) 
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case {2.4} 
m=exp (x) ; 
case {2.6,2.8} 
f£=f+exp (x<); 
case {3} 
=m+exp (x) ; 
otherwise 


£=f+exp (x) ; 


syms t 
Am exact=int (exp (tt), 1.8, 3.4) 
pretty (Am_ exact) 
error=Am exact—Am approx 
pretty (error) 
E_ t= (error/ (Am approx+error) ) *100 
pretty (E_tÈ) 
E_a= ( (Am approx—Am_ exact) /Am_ approx) *100 
pretty (E_a) 
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4-Lagrange Interpolating Polynomial Method 


Lagrange’s interpolation method uses the formula 


(X - K(X - X,)...(&X — X,) (EDE IR RK) 
i FS E TE ET ET 7 1 TET TOE TT ETN 
(Xo > K1 J(&) 7 >) ۰..(X ¬ (KX = K))J(X = X>) ...(X ¬ Xg) 
(X-XE-—Xr}...(X—XK. 1| 
۴ 0 1 r (x) 


(KX = KgJ(X, = 32) ---(X-—X,-1) 


EXAMPLE 
Given the data points 


use Lagranges method to determine y al xX = 1. 


Solution 


x) — س 2)1 ( (ڃ¥±‎ 3( 
X2) (X1 = x3) | 2( (0 > 3) 


X11) — x¥3[( 1 0J))1 —- 3( 


X11 — X2) _ Uil — Z2) 
- xa x) 3-03 2 


10= ) )x-×2( * )>×->×>3( ( / ) )×1-×2( * )×1-<×3( ( 
11= ) )×x-×1( * )>x-×>3( ( / ) )×2-×1( * )>2-×3( ( 
12= ) )×-×1( * )>×-×2( ( / ) )×3-×1( * )×3-×2( ( 
y=y0*1O+y1*1l1+y2*12 


5 کے ج سے س س ب م س س ج ج س س ج ج س س ب س سے ب سے ب م س س و س د حح س م س ص س س س س س م س کے س س ج بس سے س‎ e e e e 
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| Example 2 | 


Construct the polynomial interpolating the data by using 
Lagrange polynomials 


n SS 
E ES ٣ 
| Solution | 


ا ا ا ا ا ا ا ا ا ا ا اا اا اا ا ا ا ا ا ا ا ا ا ا ا اا ا ا ا اا ا ا اا ا ن ا ت ت ت س 


10= ) )×-×2( * )>×->×>3( ( / ) )×1-×2( * )>×1-×3( ( 
11= ) )×x-×1( * )>x-×3( ( / ) )×2-×1( * )>2-×3( ( 
12= ) )x×-×1( * )>×-×2( ( / ) )×3-×1( * )×3-×2( ( 
y=yO0*1O+y1*l1l1l+y2*12; 

collect (y) 


5 کک س س م س سے م سے س کے س ج س ص ی سے سے ب کے س ص س ص س سے س کے س ج ج ص ا ج ب کے سے م س کے س ج سے س سی س‎ e e س‎ e o n 


clc 

syms x 

p=0; 

s=]1 1/2 3]; 
E£=[3 —-10 2]; 
n=length (s) ; 
for i=1:n; 


1=1; 

for j=1:n; 
1٤ )~=[(; 

1=((x-s(]J))/(s(i)-s(j)))*l; 
end 
end 
p=l.*f (i) +p; 
end 


p=collect (Pp) 


و س س کے س ج ج ہے س سے س سے س کے س کے س ص سے س ہے س سے ب جک س سے سے کے سے ہے س کے س کے ی کے سے ہے سے کے سے کے س کے س سے کے سے کے س 5 


م. احمد محمد الفلاح الرابطى Ahmad_engineer21 @yahoo.com‏ 


الى كل من استفاد من هدا الجهد لا نسألكم الشكر ولا الثناء انما نسالكم دعوة صادقة في جوف الليل عسا ان تنفعنا في يوم تزل فيه الأقدام 


| Example 2 | 


Construct the polynomial interpolating the data by using 


Lagrange polynomials 


TTT e 
E ES ٣ 
n. 
| Solution | 


ie e e e i iE Lagrange's interpolation method-——————- 
clc 

x=input (' enter value of x:') 

p=0; 

s=]1 1/2 3]; 

£=[3 —-10 2]; 


n=length (s) ; 
for i=1:n; 


1=1; 

for j=1:n; 
ے٤‎ )تہ٣=[(;‎ 

1=))x-s)(j))/(s(i)-s(j)))*1; 
end 
end 
p=l1l.*f (i1) +p; 
end 


P; 
fprintf ('\n p(%3.3f)=%5.4f',x,P) 


س س لے ست سے سے سے سے سے سے سے کے سے سے کے کے سے سے کے کے کے کے سے سے سے سے سے کے کے سے سے کے کے سے سے کے کے سے سے سے کے سے کے کے سے کے سے کے کے کس 5 


syms x 
p=0; 
for i=1:n; 
1=1; 
for j=1:n; 
1£ )=((; 
1=))x-s(j))/(s(i)-s(j)))*1; 
end 
end 
p=l1l.*f (i1) tp; 
end 


p=collect (Pp) 


کے س کے س کے ی کے س کے سے ج س ج س ج س س س ہے جک س سے سے کے س س س ہے سے کے س کے سے کے سے کے سے سے سے کے س کے سے کے سے کے سے سے کے سے 5 


p = -283/10 -53/5 *×^2 + 419/10 *× 
enter value of x:5 
xX =5 


p)5.000(=-0 


ی س ی س نے ی س نے سے سے نے ا سے ی نے س کے ا لے نے سے ا نے نے ن ا ن کے سے کے ت س نے ا سے سے سے کے 
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| Example 3 


| Find the area by lagrange polynomial using 3 nodes | 
Solution | 


FO=6.04964; 

F1=13.464; 

F2=29. 964; 

10= ) )×-×2( * )×>×-<×3( ( / ) )×1-×2( * )×1-<×3( ( 
AO=int (10,1.8, 3.4) 

11= ) )x×-×1( * )>×->×>3( ( / ) )×2-×1( * )×>2-<×3( ( 
Al1=int (l11,1.8, 3.4) 

12= ) )×-×1( * )>×-×2( ( / ) )×3-×1( * )×3-×2( ( 
A2=int (12,1.8, 3.4) 


ا ا ا ا ا ا ا اا اا ا ا اا ا ااا ا ا ا ا ا ا ا اا ا ا ا اا ا ا ا ا ا ت ت ت ت ت س 


F=FO*AO+F1*A1+F2*AZ2 
collect (F) 


سس سے سے سے سے ہے سے سے سے سے س سے سے سے ہے سے سے سے سے سے سے سے سے سے کے کے سے س سے کے سے سے سے کے سے ہے سے سے سے سے سے کے سے کے کے کے سے سے سے سے سے سس 5 


E i e e o Lagrange's interpolation method-—-— 


syms x 

format long 

p=0; 

s=[]1.8 2.6 3.4]; 

f= [6.04964 13.464 29.964]; 
n=length (s) ; 

for i=1:n; 


1=1; 
for j=1:n; 
ے٤‎ )تہ٣=[(;‎ 
1=))x-s(j))/(s(i)-s(j)))*1; 
end 
end 


A=int (l,s(1),s (n) ) 
pP=A*f (i) +p; 
end 
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5-Mid Point Rule‏ 


Example 


Find the mid point approximation for 


using n=O 


Solution 


| Jf (x)dx= Ca^2+Ddx 
a -1 


clc 


n=6; 

h= (b-a) /n; 

£=0; 

for i=1:n; 
c=a+ (i-1/2) *h; 
£= (C21) 


6- Taylor series 


A function f(x) which possesses all derivatives up to order n at a point x = x, can be expanded 
in a Taylor series a5 


iE) 


ا 


T(x) = f(xy) + F(xX)g)J(& — xo) + (X-xXg) +... 


If x, = 0, reduces to 


"(0) 2 


f(x) = f(O)+PF(O)x + XxX + 


ل 
ك 
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Example 


Compute the first three terms of the Iaylor series expansion lor the function 


y = T(x) = tanx 


al a = r/4. 
Solution: 
The Taylor series expansion about point a is given bY 


"(a 3 


a = ۴)۵( + ۴ لع + (ھ - ٭)(4)‎ (x - a)” + (xa) +... 


and since we are asked to compute the first three terms, we must find the first and second deriva- 
tives of f(x) = tanx. 


From math tables,  tanx = sec x, SO f (x) = sec x. To find f'(x) we need to find the first 
E 
2 J : dû 
derivative of sec x, so we let z = sec x. Then, using q, SêCX = secx - tanx , We get 


dz 


d 2 
= Dsecx—sSEeECcX = ZsecxX(SecX - tanNX) = 2Zsec XK - tan&K 
dx dx 


Next, using the trigonometric identity 


1 1 
Sec X = tah xX + 1 


and by substitution , we get, 


dz ا‎ = 2(tan x + 1)Jtanx 
dx . 
Now, at point a = r/4 we have: 
| = | 1 re | ا 7 سا 1 = ا‎ = r) 2 1" e: | 1 1 کڪ ا ڪ‎ 
f (a) او‎ tanl 7 1 f(a) (7) 1+1 =2 f"(a) a 2(1 +1)1 =4 
and by substitution into (6.125), 
f(x) = 1+2] x — 1| +2x- E 
1 ) 4 1 4 | 


We can also obtain a Taylor series expansion with the MAT LAB taylor(f,n,a) function where f 
is a symbolic expression, n produces the first n terms in the series, and a defines the Taylor 
approximation about point a. 
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The following MATLAB script computes the first 8 terms of the Taylor series 
expansion of y = f(x) = tanx about a = r/4. 


Taylor s6F163 ===‏ ت کی 
clc‏ 

a=pi/4; 

sym x 

y=tan (x) ; 


z=taylor (y, 8,a); 
pretty (Zz) 


ھے ھے ھے ‏ ےد کے سے ھے سے ےھ سے سے کے سے کے سے کے سے سے سے سے کے ےھ سے سے سے سے سے سے سے کے سے کے کے پھد ےھ کے سے سد س د سے سے سے سے ےھ سے سے کد د سے سے o‏ 


Example 


Express the function 


in a Maclaurin’s series. 


Solution: 


A۸ Maclaurin’s series has the form. that is, 


I(x) = f(0) + f(0)x + x +... + 


For this function, we have f(t) = e and thus f(0) = 1. Since all derivatives are e’, then, 
f(0) = f"'(0) = f"(0) = ... = 1 and therefore, 


2 3 
E i, 
E TT 


ا 


MATLAB displays the same result. 


a E E E E E E E EE Si EE Taylor series —=————————————- 
Cl 

syms t 

fn=taylor (exp (t) ) ; 

pretty (fn) 


سے سے ہے سے کے کہ سے سے سے سے مک کے سے س س مک کک کے سے سے مک سے سے سے سے سے کے کے م جک سے سک سے مس سے س سے سے کے سے کے س س سے سے سے کے کے سے سک کے کے o‏ 
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Numerical Differentiation 
1-Finite Difference Approximations 
The derivation of the finite difference approximations for the derivatives of f (x) 


are based on forward and backward Taylor series expansions of f (x) about X, 
such as 


O O TOGO TT 
flx+MmM = fN + hf E TE ا‎ (A F< 


N 14 |‏ _ ا ) _ ا 
flx— fh) = f(x) — hf (xX) + a 1 3 (xX) + I‏ 


2) 2 


flix +2 = f(x) +2hf'(x) + ا‎ (x) + اک ا‎ ۴ 


١ (21) 7 (2R) , ع‎ 
f(x - 2R) = f(x) -2hf\( E a 


We also record the sums and differences of the series: 
: : . - 7 رر‎ " f" 4 ۳ 
flx+ M+ f(x — N =2 f(x) +R f(x) + 3 I (X) 4 ° 


e FF o. 
f(x+ mM — f(x — Rh) = 2hf (x) + f" (٭)‎ 


AM on 
flix +2 + f(x - 2R) = 2 f(x) + 4R f"(x%) + f (٭)‎ f. 


TS 
flx +28) — f(x - 2M) = 4hf (x) 1 =— f(x) 1° 


First Central Difference e 


The solution of Eq. (f) for f(4 


(+m -— fxm FH 
f] = fash f _ | (x) 


keeping only the first term on the right-hand side, we have 


flx+h— f(x — RA) 


O(h*' 
7 + Oh“) 


f) = 
which is called the first central difference approximation for f'(x). The term O(h*) 
reminds us that the truncation error behaves as j. 

From Eq. (e) we obtain 


wa _ f+ —-2 f + f(x) RF o, 
f) = 72 . f x + 
Or 
I flx+h) -2 f(x) + fx — + O(n) 


1َ 
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central difference approximations for other derivatives can be obtained from 
Eqs. (a)-h) in a similar manner. For example, eliminating /#'(x) from Eqs. (f) and (h) 
and solving for f”"(x) yield 
f"() = Fix + 2h) = 2 f(x + ا‎ mM — f(x — 2R) + O(n) 
ے2‎ 
The approximation 


fap EE re TE SPM IE oq) 


First Noncentral Finite Difference Approximations 


These expressions are called forward and backward finite difference 
approximations. 


Noncentral finite differences can also be obtained from Eqs. (a)—h). Solving 
Eq. (a) for f’(x) we get 


fFlx+h— F(X) _R 


f )x( = 7 


Keeping only the first term on the right-hand side leads to the ftrst forward difference 
approximation 


9 0 


Note that the truncation error is now O(h), which is not as good as the O(f*) error in 
central difference approximations. 


We can derive the approximations for higher derivatives in the same manner. For 
example, Eqs. (a) and (c) vield 


wa __ JAX FAM -2 f(x +R) + f) 
= h2 


+ O(B) 
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Second Noncentral Finite Difference Approximations 


Finite difference approximations of O(h) are not popular due to reasons that will be 
explained shortly. The common practice is to use expressions of O(f*). To obtain 
noncentral difference formulas of this order, we have to retain more terms in the 
Taylor series. As an illustration, we will derive the expression for] f(x). We start with 
Eqs. (a) and (c), which are 


. FF, e. aE 
f(r +B) = f(x) + hf’ (x) + > f" + = f"( + rf) +° 


. 2 TS 
flx + 2h) = f(x) + 2hF'(x) + 2h“ f" (x) + | (x) + al (xX) H+ <° 


We eliminate f“(x) by multiplying the first equation by 4 and subtracting it from the 
second equation. The result is 


fix+ 2n —=4f(x + Mh) = -3 F(X) — 2N) + FON) 4 ° 


“ے2 
a‏ 


Therefore, 


. = f{(x+2R+4f(x+R—3f(x) FH. 
f(0) = flix + +h) -3f( E E 


0 — f(x+a2amM +4 f(x +h) —3 f(X . 
a FIX + ر ا ا‎ 


Thus Equation is called the second forward finite difference approximation. 
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EXAMPLE | 
Use forward difference approximations of oh to estimate the first 
o derivative of 
fx = -0.1.*x.^4-0.15.*x.^3-0.5.*x.^2-0.25.*×+1.2 


| solution | 


% Use forward difference approximations to estimate the first 
% derivative of fx=-0.1.*x.^4-0.15.*x.^3-0.5.*x.^2-0.25.*x+1.2 
Cl 

h=0.5; 

x=0.5; 

x1=x+h 

Exx=[=0:1 =0.15 =0.5 =-0.25 1:27 

fx=polyval (fxx, x) 

fxl1=polyval (fxx, x1) 

tr va=polyval (polyder (fxx) , 0.5) 

fda= (fx1-fx) /h 

et= (tr val-—fda) / (tr val) *100 


س سے سے سے سے سے کے کے سے سے سے سے سے سے سے سے کے سے سے سے سے سے سے سے سے کے کے سے سے سے سے سے سے سے سے سے سے کے سے سے سے سے س سے سے کے سے سے سے سے سے سے سے سے سے کے کے سے ست س سے سے سے کس 6 
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| EXAMPLE 
Comparison of numerical derivative for backward difference and central 
difference method with true derivative and with standard deviation of 0.025 
x = [0:pi/50: pi]; 
yn = sin(x)+0.025 
| True derivative=td=cos(x) 


| 
| solution | 


ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ت ا ت ت ت س 9 


%< Comparison of numerical derivative algorithms 
x = [0:pi/50:pi]; 

n = length (x); 

% Sine signal with Gaussian random error 

yn = sin (x)+0.025*randn (1,n); 

% Derivative of noiseless sine signal 

td = cos(x); 

% Backward difference estimate noisy sine signal 
dynb = diff (yn) ./diff (x); 

subplot (2,1,1) 

plot (x(2:in),td(2:n),x(2:n),dynb, 'o') 

xlabel ('x') 

ylabel ('Derivative') 

axis ([0 pi -2 2]) 

legend ('True derivative', 'Backward difference') 
% Central difference 

dync = (yn(3:n)-yn(1:n-2))./(x(3:n)-x(1:n-2)); 
subplot (2,1,2) 

plot (x(2:n-1),td(2:n-1),x(2:n-1) ,dync, 'o') 
xlabel ('x') 

ylabel ('Derivative') 

axis ([0 pi -2 2]) 

legend ('True derivative', 'Central difference') 


Figure 1 
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Figure. Comparison of backward difference and central difference methods 
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Example‏ | 


Consider a Divided Difference table for points following 


0.0000 1.1487 2.7183 4.9811 


| Solution 


lI Sk! ر‎ IA Fe: Ia ا‎ Ck: ۰ Tk+3 


0.0 | 0.0000 


0.5 | 1.1487 0.8418 
3.192 0.303006 


1.5 | 4.9811 


N < IME tof o, t1 ECF fo, 1 12 EMT PIE r2) fo, 1,35, 13| 
= 000+ (r - 0.02.2074 + (r - 0.0) (r — 0.5)0.8418 + (zr - 0.0) (r - 0.5) (r — 1.00.3006 
= 2.058032 + 0.297214“ + 0.36306 
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disp (' ******** divided difference table *********') 
x=[]2 4 6 8 10] 
y=[]4.077 11.084 30.128 81.897 222. 62[ 
£00=y (1); 
for i1i=1:4 
£1(i)=(y(i+1)-y (i))/ (x(i+1)-x(i))7; 
£O01=f1 (1); 
end 
f1=[f1(1) f1(2) f1(3) f£1(4)] 
for i=1:3 
£2(i)=(f£1(i+1)-£1(i))/ (x(i+2)-x(i)); 
£O02=f£2 (1); 
end 
£2=[f2(1) £2(2) £2 (3) ] 
for i=1:2 
£3 (i)=(f£2 (i+1)-£2 (i) ) / (x(i+3) -x(i)) ; 
£03=f£3 (1); 
end 
£3=[f£3 (1) £3 (2) ] 
disp (1 KX KKK KKK KC DC KK | O DC DK KC KK KC KK KK ok o K 1) 
y=input ('enter value of y:') 
pP4x=f£00+ ( (y-x (1)) *f£01)+ ( (y-x (1)) * (y-x (2)) *£O02+ ( (y¬—x (1) ) * (y¬— 
x(2))*f£02)) 
fprintf ('\np4 (%3.3f)=%5.4f',y,P4x) 
syms Yy 
p4x=£00+ ( (y-x (1) ) *£01) + ( (y-x (1) ) * (y-x (2) ) *£02+ ( (yx (1) ) * (y- 
x(2))*f£02)) 


س ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت ت کت س کت کت کت کت کے کت کت کے کت کے کے کے کت کے کت کت کک کت کک کک کے 0 


£1 = 3.505 9:0 25. 5 70. 35 
£2 = 1. 5046 4.0906 LLL, 
£3 = 0.4310 1.1714 


p4x = -293/100+7007/2000*y+12037/4000* (y-2) * (y-4) 
enter value of y:8 

y = 8 

p4 (8.000)(=97.3200 
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| Example { H.W } | 
Find the divided differences (newten's Interpolating) for the data 


and compare with lagrange interpolating. 


| 3 | 1# * 
#* إ * | 7ة | *_ 
| مصەفutصق‏ ا 


Solution 


Xkkkkkkkkktkkkk divided difference table XXXXXXXXXXXKKK Kok 
£1 = 


26.000000000000000 4.800000000000000 


£2 = 


—-10. 600000000000000 


enter value of y:5 


p4 (5.000) =-83. 8000 
px = -283/10-53/5*y^2+419/10*y 


enter value of x:5 


Pp (5. 000)(=-83. 8000 
p = -283/10-53/5*m^2+419/10*m 
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i es ome o e ee O o Gi SE Solve IW‏ 


CLE 
disp ('******** divided difference table ****x*****') 
x=[1 0.5 3]; 
y=[]3 -10 2]; 
£00=y (1); 
for i=1:2; 
£1(i)=(y(i+1)-y(i))/(x(i+1)-x(i)) ; 
£O01=f1 (1); 
end 
£1=[f£1(1) £1 (2) ] 
for i=1; 
£2(i)=(f£1 (i+1) -£1 (1)) / (x(i+2) -x* (3) ) ; 
£O02=f2 (1); 
end 
E£2=fÊ2 (1) 


newtens Interpolating }=-——————————————=—=—-—- '(‏ { =" ) 15ل 
y=input ('enter value of y:');‏ 

px=f00+ ((y-x(1))*f£01)+ ((y-x (1)) * (y-x (2)) *f£02); 

fprintf ('\npx (%3.3f)=%5.4f',y¥, PX) 

syms Yy 

px=f00+ ((y-x(1))*f£01)+ ( (yx (1)) * (y-x (2)) *f£02); 

px=collect (px) 


£ e e o OS SO EE E COMBO WIL = 

Lagrange's interpolation method-————————————————-‏ د تي 
COMPILE WIL ma '(‏ ا ) 3p‏ 
Lagranges interpolation method-————————————————- '(‏ >== ') ح156 
m=input (' enter value of x:');‏ 

p=0; 

s=]1 1/2 3]; 

f=[3 —-10 2]; 


n=length (s) ; 
for i=1:n; 


1=1; 
for j=1:n; 
1£ )=7( ر‎ 
1=((m-s(j))/(s(i)-s(j))) *l; 
end 
end 
p=l1l.*f (i) +p; 
end 
P/ 
fprintf ('\n p(%3.3f)=%5.4f',m,pP) 
syms m 
p=0; 
for i=1:n; 
1=1; 
for j=1:n; 
ة٤ ;([=~ہ)‎ 
1=((m-s(jJ))/(s(i)-s(j)))*l; 
end 
end 
p=l1.*f(i)+p; 
end 


p=collect (Pp) 


5 ج س کے کے ب س س م س‎ e e س س ص س سے س م س کے س س س کے س ج س س س سے ج ب کے س م س کی س س س سے س مک س س س س د کے ب ج ب ص س س س‎ e e e e e 
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| Example { H.W } 


| Estimate the In(3) for | 


a) Linear Interpolation. 
B) Quardratic Interpolation 
compare between a&b 


Solution 


a) Linear Interpolation. 
F1 (x) =f (x0) + ( (f (x1) —f (x0) ) / (x1¬-x0) ) * (x—-x0) 
b)Quardratic Interpolation 


Ê2 (x) =bO+b1* (x—-x0) +b2* (x—-x0) * (x—¬-<1) 


b0= f(x0) = 0.693147180559945; 
b1= (f(x1)-f(x0))/(x1-x0) = 0.34657359027973 
b2= ((f(x2)-f(x1))/(x2-x1))-b1/(x2-x0) = -0.03596025905473; 


fxl1l = 0.693147180559945-0.346573590279973 )×-2( 
inter value x:3 
£fxl1l = 1.039720770839918 


fx2 = 0.346573590279973X+(-0.035960259056473X+0.071920518112945(*)X-4( 
inter value x:3 
£xXZ2 = 1.075681029896391 


Elt1 =5.360536964281382 ۶% 


ElLZ2Z = 2.087293124994937 >< 


Quardratic Interpolation 1s better than Linear Interpolation 
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o === Solve H.W ----------------------============ سسس‎ 
o === a) Linear Interpolation------------------------------ 
o === b) Quardratic InterpolatiOon------------------------ 
o === compare between a& D---------------------------- 
clc 


x=input('inter value x:'); 
format long 
xi=[]2 4 6]; 
fx=[log{2) 1log(4) log(6)]; 


disp(--------- a) Linear Interpolation----------------------- ( 
fx1=fx(1)+((fx(2)-fx(1))/(x1(2)-x1(1)))*(x-x1(1)) 
disp(--------- b) Quardratic Interpolation----------------- 3 
bO=fx(1); 


b1=(fx(2)-fx(1))/(x12)-xi1)); 
b2=(((fx(3)-fx(2))/(xi(3)-xi(2)))-b1)/(xi(3)-x1(1)); 
fx2=b0+b1*(x-xi(1))+b2*(x-xi(1))*(x-xI@2)); 

0 pretty(fx2) Voexpand(fx2) % collect(fx2) 


disp(---------- compare between a&b---------------------- ( 
Tv=log(3); 
disp('--------- a) Linear Interpolation----------------------- '( 
Et1=abs((T'v-fx1)/Tv)*100 
disp(--------- b) Quardratic Interpolation----------------- '( 
Et2=abs((T'v-fx2)/Tv)*100 
if Et1>Et2; 

disp('Quardratic Interpolation is better than Linear Interpolation’) 
else 

disp(' Linear Interpolation is better than Quardratic Interpolation’) 
end 
SyIHIS X 
diSp(========- a) Linear Interpolation----------------------- 3 
fx1=fx(1)+((fx(2)-fx(1))/(x1(2)-x1(1)))*(x-x1(1)) 
disp(-===----- b) Quardratic Interpolation----------------- ( 
bO=fx(1); 


b1=(fx(2)-fx(1))/(x1(2)-x1(1)); 
b2=(((f<(3)-fx(2))/(x1(3)-x1(2)))-b1)/(x1(3)-x1(1)); 
fx2=b0+b1*(x-x1(1))+b2*(x-x1(1))*(x-x1(2)) 
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The Bisection Method for Root Approximation 


we can compute the midpoint x,, of the interval x, < x < x, with 


RR 
XK mm 
7و د‎ = . 


آا. 


Knowing x, , we can find f(x, ). Then, the following decisions are made: 


1. If f(x 


This indicates that x, and x, are on the left side of the x-axis crossing as shown in Figure. 


and f(x, ) have the same sign, their product will be positive, that is, f(xy) ° f(x,) =0. 


0 
In this case, we replace x, with x, . 


f(x, ) and f(x,,) are 
both positive and thus 
their product 1s positive 


f(x, ) and f(x) are 
both negative and thus 
their product 1s pos1ti1ve 


“1 Xm A2 وھ 1ک‎ 0 


irure . Sketches to illustrate the bisection method when f(x, ) and f(x.) have same sien 
Figure . Sketches to ill he bisection method when f(x,) and f(xy) h 2 


2. If f(x, ) and f(x, ) have opposite signs, their product will be negative, that is, f(x, ) ‘f(x, ) <0. 
This indicates that x, and x, are on the right side of the x-axis crossing as in Figure. In 


this case, we replace x, with x, . 


f(x, ) and f(x, ) have f(x, ) and f(x, ) have 
opposite s1gns and thus opposite signs and thus 
theır product 1s negative their product is negative 


1 Xm n 


Figure . Sketches to illustrate the bisection method when f(x, ) and f(x, ) have opposite signs 


After making the appropriate substitution, the above process is repeated until the root we are 
seeking has a specified tolerance. T'o terminate the iterations, we either: 


a. specify a number of iterations 


b. specify a tolerance on the error of f(x) 
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Example | 
Use the Bisection Method with MATLAB to approximate one of 
the roots of 


yv = f(x) = 3x -2X + 6X -8 


by 

a. by specifying I16 iterations, and usıing a for end loop MATLAB 
program 

b. by specifying 0.00001 tolerance for f(x), and using a while end loop 
MATLAB program 


Solution: 


function y= funcbisect01 (x) ; 


VES FT g4 FP DS D2 FE ¥ ° 3 EO 2F gE = 8; 


© 


% We must not forget to type the semicolon at the end of the line 


above; 


© 


% otherwise our script will fill the screen with values of y 


س 


call for function under name funcbisect01.m 


disp (' xm fm')  % xm is the average of xl and x2, fm is 
1ل‎ 5p ) "= '( % insert line under xm and 


for k=1:16; 
f1=funcbisect01 (x1); f£2=funcbisect01 (x2) ; 
xm= (x1+x2) / 2; fm=funcbisect01 (xm) ; 
fprintf ('%9.6f %13.6f \n', xm, fm) % Prints xm and fm on same 
line; 

if (f1*fm<0) 

x2=xm; 
else 


xXx=1:0.05:2; 

VETER DST EGS OSO EOCGO TREE 
plot (x,y) 

grid 
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Fiaure 1 - 


File Edit Yiew Insert Tools Desktop Window Help 


GOFF kh AAS # OS AG 


1 


لا 


1 


bU 
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function y= funcbisect01 (x) ;‏ 
ya 3 FRE ° Da EE 6% 3F 6 E E = 8‏ 
We must not forget to type the semicolon at the end of the line‏ % 
above;‏ 
otherwise our script will fill the screen with values of y‏ % 


م و کے س سے ب سے س سے ج س کے س ج بس س ب کے س کے س ج س س س ج س سے ج ج ب سے ب ص س س ص س سے س کے ج ج س س س س ب س ب ج ب ص س ج س س س ج س کے س بے سے س سے س ١‏ 


call for function under name funcbisect01.m 


X2=2; 


while (abs (x1-x2) >2*tol) ; 
f£1=funcbisect01 (x1) ; 
£2=funcbisect01 (x2) ; 

xm= (x1+x2) /2; 

f£m=funcbisect01 (xm) ; 

fprintf ('%9.6f %13.6f \n', xm, fm); 
if (f£1*fm<0O) ; 


xXx2=xm; 

else 

x1=zxm; 

end 

end 

1p '( 
xm fm 

1. 500000 17.0.0 
1.250000 4.7149023 
1.12500 1. 308441 
1. 062500 0.03838 
1. 0 -0. 506944 
1.046875 -0. 4 
1.054688 -0. 15 
1.058594 -0.N 5 
1. 060547 0. 002604 
1.059570 -0.015168 
1. 060059 —-0. 006289 
1. 060303 -0.004 
1. 060425 0.000380 
1. 060364 -0.000732 
1. 060394 —-0.000176 
1. 060410 0.000102 
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Example | 
Use the Bıisection Method with MATLAB to approximate one of 
the roots of (to find the roots of) 
Y=f(x)= x.^3-10.*x.^2+5; 
That lies 1n the 1nterval ( 0.6,0.8 ) by specifying 0.00001 tolerance 
for f(x), and using a while end loop MATLAB program 
Solution: 


ا ا ا ا ا ا ا ا ا ااا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ت ت ت ت ت ت س 9 


function y= funcbisect01 (x); 
y = x.^3-10.*x.^2+5; 


© 


% We must not forget to type the semicolon at the end of the line 
above; (% otherwise our script will fill the screen with values of y) 


ا ا ا ا ا اا ا ااا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا ا اا ا ا ا ا ا ا ا ا ا ا ن ا ت ت ت س 9 


call for function under name funcbisect01.m 


clc 

x1=0.6; x2=0.8;tol=0.00001; 

disp (' ج کے‎ '( 
disp (' xm fm'); 
disp (' س س ت‎ '( 


while (abs (x1-x2) >2*tol) ; 
f£1=funcbisect01 (x1) ; 
£2=funcbisect01 (x2) ; 

xm= (x1+x2) /2; 

f£m=funcbisect01 (xm) ; 

fprintf ('%9.6f %13.6f \n', xm, fm) ; 
1f (f£1*fm<0O) ; 


X2=xm; 

else 

x1=xm; 

end 

end 

J15 ) جت‎ '( 
xm fm 
0.700000 0.440 
0.750000 --5 
0.725000 0.18 
0.737500 - 2 
0.731250 0.04733 
0.734375: 7 
0.735938 -.53 
0.735156 -0.007228 
0.734766 -0 
0.734570)... 4 
0.734668 -0.00083 
0.734619 -0.004 
0.734595: 5 
0.734607 -l - 5 
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Newton-Raphson Method 


The Newton-Raphson formula can be derived from the Taylor series expansion 
of f(x) about x: 


f (xi) = f(D + f (D(x — x) + Ox — x; (a) 
Ifx;,ر‎ 1s a root of f(x) = 0, Eq. (a) becomes 
0= f(x) + f (xo) (x1 — x) + O(%+1 — X)* (b) 


Assuming that x, is a close to X;,ı, we can drop the last term in Eq. (b) and solve for 
X;,1. The result 1s the Newton-Raphson formula 
ا‎ )( (c( 
(x) 
Ifx denotes the true value ofthe root, the error in X; iS E; = X — x;. Itcan be shown 
that If xX, 1s computed from Eq. ( c ), the corresponding error 1s 


1 = Xi 


indicating that the Newton-Raphson method converges quadratically (the error 1s the 
square of the error in the previous step). As a consequence, the number of significant 
figures Is roughly doubled in every iteration, provided that x; 1s close to the root. 


Tangent line. 


| Figure( a )Graphical interpretation of the Newton-Raphson 
۱ f(x) formula. 

| 

| 


A graphical depiction of the Newton-Raphson formula is shown in Fig. ( a ) The for- 
mula approximates f(x) by the straight line that 1s tangent to the curve at X;. ThUS Xi, 1 
is at the Intersection of the x-axis and the tangent line. 

The algorithm for the Newton-Raphson method 1s simple: it repeatedly applies 


Eq. ( c ), starting with an initial value xo, until the convergence criterion 
xy — X1| < € 


is reached, e being the error tolerance. Only the latest value of x has to be stored. Here 
is the algorithm: 


1. Let x be a guess for the root of f(x) = 0. 
2. Compute AX = — f(X/PF(). 
3. Let x > x + Ax and repeat steps 2-3 untill |AX| < €. 


م. احمد محمد الفلاح الرابطى Ahmad_engineer21 @yahoo.com‏ 


الى كل من استفاد من هدا الجهد لا نسألكم الشكر ولا الثناء انما نسالكم دعوة صادقة في جوف الليل عسا ان تنفعنا في يوم تزل فيه الأقدام 


EXAMPLE 


A root of f(x) = x3 — 10x* + 5 = 0 lies close to x = 0.7. Compute this root with the 
Newton-Raphson method. 


Solution 


The derivative of the function 1s f(x) = 3x* —- 20x, so that the Newton- 
Raphson formula in Eq. ( c ) is 


f) x —-10x% +5 2x" -10% 5 


a TT OT Î o Û GET TE RFT‏ س 
3x2 - 20x x )3x - 20(‏ )۴(0 4 


It takes only two iterations to reach five decimal place accuracy: 


20.7) - 10)0.7(* 5 
U7 EOT = 20 


= 0.735 6 


5 - 336(2 10)0.735 - )20.73536 ۰ 
ا ڪڪ ي ٣‏ 
|20 - (36 3(0.735] 36 0.735 


= 0.134 60 
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| Example | 


Use the Newton—-Raphson Method to estimate the root of f(<)=e^(-x)-x, 
employing an Initial guess of x0=0 


x= [0]; 
tol=0. 0000000007; 
format long 
for 1=1:5; 
fx=exp (=x (1) ) x (i) ; 
fxx=-exp (-x(1))-1 ; 
Efxxx=exp (-x(i1)); 
x(1+1) =x (i) - (fx/Êxx) ; 
T.V(i)= (abs ((x(i+1) -x (i) ) /x(i+1) ) ) *100; 


end 
for 1i=1:5; 
e(i)=x(6)-x(i); 
fxx=-—exp (-x (6) )—-1 ; 
fxxx=exp (-x (6)); 
e(i+1) = (-fxxx/2*fÊxx) * (e (i) ) ^2; 
end 


1f abs (x(i+1)-x (i1) ) <tol 
disp (' enough to here') 
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enough to here 


0.500000000 
0.566311003197218 
0.567143165034862 
0.567143290409781 
0.567143290409784 


1.0e+002 * 


1.000000000000000 
0.117092909766624 
0.0014675 
0.000000221063919 
0.000000005 


0.567143290409784 
0.067143290409784 
0.000832285 66 
0.000000125 7422 
0.000000003 
0.000000000 
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The secant Formula Method 


A popular method ofhand computation is the secant formula where the improved 
estimate of the root (x, , ) is obtained by linear interpolation based two previous 
estimates (x, and Xx): 


1 = ا 


F(x; )‏ ك 
Fon = FON‏ 


M+] = 


| Example 


Use the The secant Formula Method to estimate the root of 
{(x)=e^(-x)-x, employing an initial guess of x(1-1)=0 & x(0)=0 


E 


Solution 


e as The: SéCanNt FOULS MElINOM acini 
CLE 
x=[0 1]; 


eT 1 = X1 
ہے = [ لإا‎ f(x 8 Fx) 


TV=0.567143290409784; 
format long 
for 1=2:6; 
Efx=exp (-x(i1-1) )-x(i-—-1); 
fxx=exp (-x(1))-x<(i1); 
x(1+1) =x (i) —-((x(1i) -x (i-1) ) *fÊxx) / (fxx—fÊx) ; 
E T(i)=(abs ( (TV-x(i+1) ) /TV) ) *100; 
end 
I155) == '( 
disp(' X(i+1) ') 
JI1SD (mau '( 
جد‎ 
E o e e E E '( 
disp (' E 7T '( 


م. احمد محمد الفلاح الرابطى Ahmad_engineer21 @yahoo.com‏ 


الى كل من استفاد من هدا الجهد لا نسألكم الشكر ولا الثناء انما نسالكم دعوة صادقة في جوف الليل عسا ان تنفعنا في يوم تزل فيه الأقدام 


. 000000000000000 
. 610962 
. 563838314 
. 567170358415 
. 5671433066024963 
. 567143290409705 


OOO Oop 


. 032634281468 
. = 442 
.N N 440 
. 000002855570996 
. 00000000001397 


O O O O oO 
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| Example | 


Use N.R. Quadratically Method to estimate the multiple root of 
1(x<)=x^3-5x^2+7x-3, initial guess of x(0)=0 


f) f" )x( 


+1 س‎ xX = 


f(x)— fx x) 


Solution 


The N.R. Quadratically Method _ ————————————-————-‏ متتو 
Cle‏ 

TV=1; 

x= [0]; 

format long 

for 1=1:6; 


E£x=x (1) ^3-5*x (1) ^2+7*x (1) -3 
E£xx=3*x (1) ^2—-10*x (1) +7 

Efxxx=6*x (1) ~10 

x (1+1) =x (i) - (fx*fÊxx) / ((fxx) ^2—-fx*fÊxxx) ; 
E T(i)=(abs ( (TV-x(i+1) ) /TV) ) *100; 


end 

es es i SS '( 

disp(' X(i+1) ') 

e E) 2 e e E '( 

جد 

JIS (mau '( 

disp (' E T '( 

J155 ) == '( 

E T' 

Jip (1 ===> '( 
e Se E EE Multiple Roots-———————- 
>--٤×= )×-3( )×-1( )×-1( = 
clc 


for x=—-1:0.01:6; 
E£x=x.^3-5.*x. ^247. *x¬-3 
plot (x, fx) 
hold on 

end 

grid 

title (' (x—-3) (x—-1) (x—-1) ') 

xlabel ('x' ) 

ylabel (' f£x') 
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1.10520637 
1.003081664098603 
1.000002381 6 
1.44444 A A 12 
1.M0MMشھ‎ A A 5 
1.00 | ۇ$‎ 5 


1.334333 73006 
0.308166403 
0.00238149381548 
0.114445 
0.000000075 
0.000000075 


Fıgqure Î 
Fils Edl Yew Inset Tools Desklap Window Help 


TE HêR| AA S| #| OE a A 
[ Fx=t-Ax-Tix-1 ] 
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| Example | 


Use the Newton—Raphson Method to estimate the root of 
(x)=x^3-5x<^2+7x-3, Initial guess of x(0)=4 


x=[4]; 
tol=0. 0007; 
TV=3; 
format long 
for i=1:5; 
fx=x (1) ^3-5*x (1) ^2+7*x (1) -3; 
fxx=3*x (1) ^2-10*x (1) +7; 
x (1+1) =x (i) -(fx/fxx) ; 
E T(i)=(abs ( (TV-zx(i+1) ) /TV) ) *100; 
end 
for i=1:5; 
e(i)=x(6)-x(i); 
Ê£fx=x (1) ^3-5*x (1) ^2+7*x (1) -3; 
E£xx=3*x (1) ^2—-10*x (1) +7; 
E£xxx=6*x (1) ~10; 
e(i+1) = (-fxxx/2*fxx) * (e (i1) ) ^2; 
end 


1f abs (TV-x (1+1) ) <tol 
disp (' enough to here') 
Jip (awa '( 
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enough to here 


4.000000000000000 
3. 400000000000000 
3. 100000000000000 
3.008695653 
3. 0000746409192 
3. 000000005570623 


3. 0 
CEECECECCCECEKECERKE KECE 
0.289855072463781 
0.0024880263730 60 
0.000000185687436 
0.000000007462475 


-0. 999999994427 
-0. 399999994427 
-0.09999999442 7 
-0.0086956466030 
-0.000074635 9 
—-0.000000089144954 
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Gauss Elimination Method‏ 


س 


| Example | 
Use the Gauss Elimination Method with MATLAB to solve the 


followıIng equations 


2X1 +X2-X3=5--------==-==--=========== (1( 
X1+2<2+4<3=10---------=============- (2) 
EO PI e E (3) 
| Solution: 
جي‎ GaUS5S ElinNination Neltigog—=——————= ج‎ 
clc 
A=[2 1 -1;1 2 4;5 4 —-1]; 
b=[5 10 14]; 


1f size(b,2) > 1; b = b'; end % b must be column vector 
n = length (b); 

for k = 1:n-1 % Elimination phase 

for i= kt+1:n 

1f A(i,k) ~= 0 

lambda = A(i,k)/A(k,k); 

A(i,k+1:n) = A(i,k+1:n) — lambda*A(k,k+1:n); 
b(i)= b(i) ¬ lambda*b (k) ; 

end 

end 

end 

1f nargout == 2; det = prod (diag (A) ); end 
for k = n:-1:1 % Back substitution phase 
b(k) = (b(k) - A(k,k+1:n) *b(k+1:n))/A(k,k); 


fprintf (' 
end 
x = b 
X= 
4 
-1 
2 
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ترقبوا المزيد من الشروحات للأمثلة فى التحليل العددى 
والرياضيات والتحكم الألى والأتصالات وإلكترونات القدرة 
ونظم التشغيل والدارات التماثلية والنظم الرقمية واسس 
الألكترونات وغيرها من المواد فى اغلب التخصصات راجين 
والعصمة. 


وفى الختام نسأل الله التوفيق والسعادة لى ولكم في الدنيا والأخرة . 
ونسأل الله الهمة في طلب العلم وبدله . واللهم صلى على النبى المصطفى 
وال بيته الطاهرين والصحابه والتابعين وتابع التابعين ومن تبعهم 
باحسان الى يوم الدين. والسلام عليكم ورحمة الله وبركاته. 


م. احمد محمد الفلاح الرابطى Ahmad_engineer21 @yahoo.com‏ 


This document was created with Win2PDF available at http://www.win2pdf.com. 
The unregistered version of Win2PDF is for evaluation or non-commercial use only. 
This page will not be added after purchasing Win2PDF. 


